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CN , Abstract 

^ . Finite difference scliemes are here solved by means of a linear matrix equa- 

1—5 ', tion. The theoretical study of the related algebraic system is exposed, and 

' enables us to minimize the error due to a finite difference approximation, while 

' building a new DRP scheme in the same time. 
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g : 1 Introduction: Scheme classes 

^ ■ We hereafter propose a method that enables us to build a DRP scheme while mini- 

^ . mizing the error due to the finite difference approximation, by means of an equivalent 

! matrix equation. 

O ■ 
cn ; 

y—^ • Consider the transport equation: 

OO , (^u, Su 

O. — + c— = , X e [0,L], t e [0,T] (1) 

: ^ 

I with the initial condition u{x,t = 0) = uo{x). 

Proposition 1.1 A finite difference scheme for this equation can be written under 
the form: 

a u,"+7 u,"''+6 Mi+r+e Mi_i"+C Ui.i"^^+^ m+i"-' = 

(2) 

where: 

ur = u{lh,mT) (3) 

/ G {i — l,i,i + l},m G {n — l,n,n + l},j = 0,...,nx,n = 0,...,nt,h,T denoting 
respectively the mesh size and time step (L = Ux h, T = ritr). 
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The Courant-Friedrichs-Lewy number [cfl) is defined as a = cr/h . 

A numerical scheme is specified by selecting appropriate values of the coefficients 
a, P, 7, 6, e, (, r], 6 and in equation ([2]), which, for sake of usefulness, will be 
written as: 

a = a^ + at , P = Px + Pt , 1 = Ix + It , 5 = 5x + St , e = ex + et , (4) 

where the denotes a dependance upon the mesh size h, while the "t" denotes a 
dependance upon the time step r. 



The number of time steps will be denoted rit, the number of space steps, Ux- In 
general, ^ rit. 

In the following: the only dependance of the coefficients upon the time step r existing 
only in the Crank-Nicolson scheme, we will restrain our study to the specific case: 

at = ^t = C = V = d = ^ = Q (5) 

The paper is organized as follows. The building of the DRP scheme is exposed in 
section [2l The equivalent matrix equation, which enables us to minimize the error 
due to the finite difference approximation, is presented in section [3l A numerical 
example is given in section ??. 

2 The DRP scheme 

The first derivative || is approximated at the node of the spatial mesh by: 

du 

{-Q^)i^Px Ui+Z" + 6x Ui+i+i" + ex ui+i^i"" (6) 

Following the method exposed by C. Tam and J. Webb in the coefficients Px, 
Sx, and Ex are determined requiring the Fourier Transform of the finite difference 
scheme (E]) to be a close approximation of the partial derivative 
(E]) is a special case of: 

du 

( — )/ ~ u(x + ih) + 6x uix + ii + l)h) + u(x + (i — 1) h) (7) 
ox 

where x is a continuous variable, and can be recovered setting x = I h. 

Denote by u the phase. Applying the Fourier transform, referred to by ^, to both 

sides of (I7j), yields: 

juuc^{Pxe' + 6xe^'^'' + exe-^'^''} u (8) 
j denoting the complex square root of —1. 
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Comparing the two sides of (E]) enables us to identify the wavenumber A of the finite 
difference scheme ([6]) and the quantity j + 5^ e-^'^'' + e"-''^'^ }, i. e.: The 

wavenumber of the finite difference scheme (E]) is thus: 

J=-j {Pxe'^ + Sxe^'^^ + exe-^^^} (9) 

To ensure that the Fourier transform of the finite difference scheme is a good ap- 
proximation of the partial derivative ( |^ )/ over the range of waves with wavelength 
longer than 4h, the a priori unknowns coefficients jS^, 5x, and must be choosen 
so as to minimize the integrated error: 



S = J\\Xh- Xh\^d{Xh) 



The conditions that £^ is a minimum are: 

d£ 



d(3x 



dSx 



d£_ 

dsx 



and provide the following system of linear algebraic equations: 

2'Kh(3x + A{h6x + hex-l) = 
Ah/3x + TT{2 5x-l) = 
4:hp^ + 27rhex = 

which enables us to determine the required values of P^, ^x, and 



Px 


= i^r 


Sx 


= 


Sx 


— popt 
^x 



/i(7r2-8) 
1 2 
2 



/i(7r2-8) 
2 



(10) 



(11) 



:i2) 



(13) 



/i(7r2-8) 



3 The Sylvester equation 

3.1 Matricial form of the finite diff'erences problem 

Theorem 3.1 The problem ([2]) can be written under the following matricial form: 

M1U + UM2 + C{U) = Mo (14) 

where Mi and M2 are square matrices respectively — 1 by — 1, nj by n^, given 
by: 



Ml 



//3 


s 


... \ 




( 7 


... \ 


e 








a 









■-. ■•. 


M2 = 


■•. 


■•. ■•. 






■•. (3 6 








V 




e p J 




[0 ... 


a J 



(15) 



the matrix Mq being given by: 
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Mo 



- e Mq — ?7 — & - 



-7<,-2 -'?«°,-2 -''^^-l 







-■5"^^ -c^L -''"L 



(16) 



and where £ is a linear matricial operator which can be written as: 

£ = £1 + £2 + /^s + A 
where £1, £2, £3 and £4 are given by: 



A(t/) = C 



V 



Un 







/ 



/ 

u\ 




m1 



(17) 



\ 



Tlx 2 ^'i^x 2 



;i8) 



<-2 / 



/ 



£3(t/) = 9 



ul 



ul 



^Tlx 2 ^TLx 2 








■ur-o y 



/o 





C4{U) = ^ 



U2 U2 
1 2 

Uq ui 



<-i 



nt-1 
2 

nt-1 



rit — 1 



(19) 



/ 



Proposition 3.2 The second member matrix Mq bears the initial conditions, given 
for the specific value n = 0, which correspond to the initialization process when 
computing loops, and the boundary conditions, given for the specific values i = 0, 
i = Ux- 

Denote by Uexact the exact solution of ([I]). 
The corresponding matrix Uexact will be: 

Uexact [Uexacti ] l<i<nx — l,l<n<nt (^0) 

where: 



Uexacti Uexacti^-^iy ^n) (^-^) 



with Xi = i h, tn = n T. 
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Definition 3.3 We will call error matrix the matrix defined by: 

E = U- Ue^act (22) 

Consider the matrix F defined by: 

F = Ml Uexact + Uexact M2 + C{Uexact) ' Mq (23) 

Proposition 3.4 The error matrix E satisfies: 

M1E + EM2 + C{E) = F (24) 

3.2 The matrix equation 

Theorem 3.5 Minimizing the error due to the approximation induced by the nu- 
merical scheme is equivalent to minimizing the norm of the matrices E satisfying 

m- 

Note: Since the linear matricial operator C appears only in the Crank-Nicolson 
scheme, we will restrain our study to the case C = 0. The generalization to the case 
£ 7^ can be easily deduced. 

Proposition 3.6 The problem is then the determination of the minimum norm 
solution of: 

MiE + EM2 = F (25) 
which is a specific form of the Sylvester equation: 

AX + XB = C (26) 

where A and B are respectively m by m and n by n matrices, C and X, m by n 
matrices. 



3.3 Minimization of the error 



Calculation yields: 

' M^M H^nn(f + {S + s) \ f + P {S + s) \ 



(27) 



The singular values of Mi are the singular values of the block matrix ( 



1. e. 



I (2/?2 + s^ + e^-{6 + e) ^^413^ + 6^ + - 26 e) 



(3^ + 6^ (3 {6 + e) 
(28) 
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of order ^^^-'-j and 



1 (2/?2 + ^2 ^ ^2 ^ ^ ^) V4/52 + 52 + e2_25£) (29) 



of order 



The singular values of M2 are a^, of order ^, and 7^, of order y- 
Consider the singular value decomposition of the matrices Mi and M2: 



M2 




(30) 



where Ui, Vi, U2, V2, are orthogonal matrices. Mi, M2 are diagonal matrices, the 
diagonal terms of which are respectively the nonzero eigenvalues of the symmetric 
matrices Mi^Mi, Ma'^Mg. 

Multiplying respectively [25] on the left side by '^Ui, on the right side by V2, yields: 

MiEV2 + U^ EM2V2 = U^ FV2 (31) 
which can also be taken as: 

^f/i Ml Vi ^Vi E V2 +^ Ui E ^U2 ^U2 M2 V2 = F V2 (32) 

Set: 

'V,EV2={ |ii ^]jU,E^U2=( ^ ^ I (33) 

^21 F22 J y E21 E22 

U,FV2={ p p 1 (34) 

1*21 r22 

We have thus: 



Ml En MiEu \ _^ ( 1 _ ( ^1 ^2 

j 1 ^^^2 I V ^21 F22 



(35) 



It yields: 



MiEi^+EiiMa = En 

MiEu = Fi2 (36) 



E21M2 = Fs 



21 



One easily deduces: 



£"12 = Ml ^ F12 
E21 = F21 M2 



(37) 



The problem is then the determination of the £^11 and En satisfying: 
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Ml Eu + Eu M2 = Eu (38) 

Denote respectively by eTj, efj the components of the matrices E, E. 
The problem [381 uncouples into the independent problems: 
minimize 

Gij + Cij (39) 

under the constraint ^ ^ 

Mui ^ + ^ = Euij (40) 

This latter problem has the solution: 



— Miii Fiiij 

Cji — ■ 2 2 



(41) 



The minimum norm solution of [23 will then be obtained when the norm of the 
matrix En is minimum. 

In the following, the euclidean norm will be considered. 
Due to ([3] 



lli^iill < lli^ll < \\Ui\\ \\E\\ IIV2II < ||f/i|| ||V^2|| \\M,U,^act + U,^actM2-Mo\\ (42) 

Ui and V2 being orthogonal matrices, respectively — 1 by ra^. — 1, rit by rit, we 
have: 

||f/if = n,-l , \\V2f = nt (43) 

Also: 

\\M,r = !^l^(2P' + S' + e') , ||M2||2 = |(a2 + f) (44) 



The norm of Mq is obtained thanks to relation ([T6 
This results in: 



llFiill < Vnt {n,-l) |||f/e..acdl y^V^^^+T^) + \\Mo 

_ (45) 
can be minimized through the minimization of the second factor of the right- 
side member of fj^3|) . which is function of the scheme parameters. 



||t^exarf|| is a constant. The quantities y a/2 + 5^ + e^, \/ a'^ + •y'^ and ||Mo 

being strictly positive, minimizing the second factor of the right-side member of (j45 
can be obtained through the minimization of the following functions: 



/2(«,7) = (46) 
/3(a,/?,7,5,e) = II Mo II 
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I.e.: 



Setting: 



/2(a,7) = V^^l + % 
f3{a,P,7,S,e) = \\Mo 



13. 


= 




= 




— popt 

^x 



(47) 



(4J 



one obtains the DRP scheme with the minimal error through the minimization of: 



gi{pt,5t,et) 
g2{ax,lx) 



0^1 +lx 

Moll 



(49) 



4 Numerical results 

We denote by t the non-dimensional time parameter. Figured] displays the norm 
of the error for an optimized scheme (in black), where (3^, ^x, ^x are given by (fT3ll . 
and a non-optimized one: numerical results perfectly fit the theoretical ones. 
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Figure 1 : norm of the error for the optimized scheme (in black) and the non optimized 
one (dashed curve). 

Figure [2] displays the L°° norm of the error for the above optimized scheme (in 
black), a seven-point stencil DRP scheme (in gray), and the FCTS scheme (dashed 
plot). As time increases, the optimized scheme yields, as expected, better results 
than the FCTS one. Also, for 15 < ^ < 100, results appear to be better than those 
of the classical DRP scheme. For large values of the time parameter, both latter 
schemes yield the same results. 

Figure [3] displays the norm of the error for the above optimized scheme (in black), 
the seventh-order DRP scheme (in gray), and the FCTS scheme (dashed plot). As 
expected, results coincide. 
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Figure 2: L°° norm of the error for the optimized scheme (in black), the DRP scheme (in 
gray), and the FCTS scheme (dashed plot). 
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Figure 3: norm of the error for the optimized case (in black) and the optimized case 
(dashed curve). 

5 Conclusion 

The above results open new ways for the building of DRP schemes. It seems that 
the research on this problem has not been performed before as far as our knowledge 
goes. In the near future, we are going to extend the techniques described herein to 
nonlinear schemes, in conjunction with other innovative methods as the Lie group 
theory. 
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